NOAA Nearshore Fish Atlas of Alaska

Data Wrangle and Exploratory Analyses, Part II

Built with R version 4.4.2.

1 Introduction

The purpose of this document is to make open and share-able the research methods used for my dissertation on nearshore fish communities in Alaska working towards a PhD in marine biology at the University of Alaska Fairbanks. Here, I cover second-stage steps of the data preparation for my second and third chapters concerning spatial and temporal distributions of nearshore fishes across the state. This and other files can be accessed via the Kachemak Bay National Estuarine Research Reserve’s github repo.

In this document I share exploratory data analyses conducted on the NOAA Nearshore Fish Atlas (NFA) database. The NFA database can be found here, https://alaskafisheries.noaa.gov/mapping/sz/. Namely, this document contains the steps taken after cleaning/wrangling raw data. In particular I produce various visualizations of the data in space and time. Part I of these data preparation steps can be accessed at https://rpubs.com/chguo1/1188614, and the resulting data objects are loaded below (file “nfaa_1.rda”).

1.1 Set up

Load required packages, define directory, set options, source/load files:

# Packages
library(tidyverse)
library(lubridate)
library(here)
library(leaflet)
library(RColorBrewer)
library(vegan)
library(ggsignif)
library(ggpubr)

# Directory
wd = here()
dirs = wd %>% list.files() %>% str_subset(pattern = "^README|^LICENSE|.md$|.Rproj$", negate = TRUE)
for (i in seq_along(dirs)) {
  name = str_replace_all(dirs[i], "^", "dir.")
  path = str_replace_all(dirs[i], "^", str_c(wd, "/"))
  assign(name, path)
  rm(name, path, i)
}

# Options

# Source/Load
load(file.path(dir.data, "nfaa_1.rda")) # wrangled data

1.2 Background and objectives

Based on my own research and related literature, I can broadly interpret these beach seine data in a spatiotemporal context. Depending on when (time of year) and where (location and habitat), I have a general structure in mind of the diversity of the community. Depending on who’s there, I have a rough guess of the relative abundance of each community member (at least for the more commonly caught ones). My hope with working with the NFA data is that this broad understanding can be formally tested with inference and statistics, and that more informative research questions can be answered.

Findings from other nearshore fish researchers largely agree that seasonality (or some related environmental condition) is a strong predictor of community or species presence. Interannual differences can exhibit wide variability and should be accounted for whenever possible. And at larger time-scales (e.g., decadal), communities exhibit changes which may relate to regime shifts occurring more broadly than the nearshore.

Spatially, large-scale (regional) effects appear to be as or more important than sub-regional or local-scale effects. However, many studies also find significant relationships in community response at these smaller scales. Likely, the most appropriate spatial consideration depends on the question being asked. So maybe a more appropriate question to ask of the NFA data is if there any significant spatial scales discernible in the data, and also how should we address those scales in future research, such as studies on subsets of the taxa or habitat management considerations.

I’ll start by exploring the structure of our ‘visits’ dataframe, and then I’ll move on to visualizing the information derived from out ‘catch’ dataframe.

2 Exploring visit data

2.1 Summary of visits

Remember that in the data wrangling stage, we created a new sample identifier based on site and date which we called VisitID. Let’s take a look at the variables associated with each VisitID:

glimpse(visits)

## Rows: 1,867
## Columns: 7
## $ VisitID    <chr> "1_2001-07-24", "1_2002-03-25", "1_2001-04-26", "1_2002-07-…
## $ Replicates <int> 5, 7, 4, 4, 5, 3, 4, 2, 6, 1, 2, 2, 6, 4, 4, 4, 4, 3, 2, 2,…
## $ Date       <date> 2001-07-24, 2002-03-25, 2001-04-26, 2002-07-24, 2003-03-22…
## $ Lat        <dbl> 58.56626, 58.56296, 58.56408, 58.56755, 58.55804, 58.56858,…
## $ Lon        <dbl> -134.9105, -134.9082, -134.9137, -134.9133, -134.9102, -134…
## $ MeshSize   <dbl> 3.2, 3.2, 3.2, 3.2, 3.2, 3.2, 3.2, 3.2, 3.2, 3.2, 3.2, 3.2,…
## $ Cluster    <fct> 327, 327, 327, 327, 327, 327, 327, 327, 327, 327, 327, 326,…

During the wrangling stage, I decided that Replicates and MeshSize should be the main covariates that needed consideration. At this point they are classed as integers/numbers, but we will probably want them to be factors instead. Cluster is a byproduct of clustering some sites together, so it may not be needed if we’re using Lat/Lon as our primary explanatory variable. Still, it could be a useful grouping factor later on. Date will likely be mutated into multiple additional periods (day of year, month, year) to see if any of those factors are useful in explaining fish data. We’ll address edits to these variables as they come up in our exploration.

2.2 Spatial distribution

leaflet(visits) %>% 
  addTiles(options = tileOptions(minZoom = 3.5,
                                 zIndex = 0.5)) %>%
  setView(lat = 65, lng = -152, zoom = 3.5) %>%
  addSimpleGraticule(interval = 5) %>%  
  addCircleMarkers(lng = ~Lon, lat = ~Lat,
                   stroke = FALSE,
                   popup = ~VisitID,
                   clusterOptions = markerClusterOptions())

In the map above we can view all of our samples from across Alaska. Leaflet options are nice because you can optionally group samples on the map and show a boundary polygon of the area they cover when you hover over their icon. At the minimum zoom level (3.5), we see the largest group spanning nearly all of Southeast Alaska, followed by a group in Southcentral that covers Prince William Sound, Cook Inlet, and the northern half of Kodiak Island. The third largest group forms a much smaller polygon area around Utqiagvik. From there, we see even smaller groups from the Eastern Beaufort Sea around Kaktovik, the Bering Strait and Southern Chukchi Sea near Kotzebue, and Bristol Bay. There are also two groups from the Alaska Peninsula and Aleutian Islands, one focused around Adak and another spanning from Sand Point to Unalaska. We can see how subsets of the samples change by zooming in. When there are single samples in view, hovering over those will tell us its VisitID name.

Obviously, this visualization of the data is not perfect and pretty rough- note that the groups are formed by the overlap of samples based on pixel radius (10). For example, the polygon of samples from Bristol Bay at the minimum zoom level actually crosses the Alaska Peninsula to include a sample from Aghiyuk Island. Although, the map is useful in gaining a general understanding of the existing clusters of samples within the NFA database.

2.3 Regions

We may want to define a spatial factor based on these rough groupings as a starting place to see if the smaller, more isolated groups should be lumped with other groups or not. Actually, the NFA already has a Region classifier per site, but the reason I do not want to use it is because they linked Region to SiteID. When we created our new VisitID, some samples now contained two Region labels- these were mostly cases from around Utqiagvik where some were labelled as Chuckchi and others as Beaufort but actually occurred on the same day and within a very short distance of each other.

Instead of using the the given regional classes, let’s just make our own variable called ‘Region’ and add it to the our visits df:

# Start with events so that we can make use of Location info for unique cases,
regions = events %>%
  select(EventID, VisitID, Lat, Lon, Location) %>%
  mutate(Region = case_when(Lat > 69 & Lat < 80 & Lon > -150 & Lon < -140 ~ "Beaufort East",
                            Lat > 69 & Lat < 80 & Lon > -160 & Lon < -150 ~ "Chukchi/Beaufort",
                            Lat > 65 & Lat < 69 & Lon > -170 & Lon < -160 ~ "Chukchi South",
                            Lat > 57 & Lat < 60 & Lon > -165 & Lon < -155 ~ "Bristol Bay",
                            Lat > 50 & Lat < 55 & Lon > -180 & Lon < -175 ~ "Aleutians Adak",
                            Lat > 50 & Lat < 55 & Lon > -170 & Lon < -165 ~ "Aleutians Unalaska",
                            Lat > 55 & Lat < 65 & Lon > -155 & Lon < -145 ~ "GOA Southcentral",
                            Lat > 53 & Lat < 61 & Lon > -148 & Lon < -130 ~ "GOA Southeast"))

# Check which samples we missed,
filter(regions, is.na(Region))

## # A tibble: 7 × 6
##   EventID VisitID           Lat   Lon Location                            Region
##     <dbl> <chr>           <dbl> <dbl> <chr>                               <chr> 
## 1    6635 1564_2006-06-22  52.7 -171. Yunaska Island                      <NA>  
## 2    6636 1565_2006-07-26  56.2 -157. Semidi Islands, Aghiyuk Island      <NA>  
## 3    6637 1566_2006-07-31  55.1 -160. Shumagin Islands, Bendel Island     <NA>  
## 4    6638 1566_2006-07-31  55.1 -160. Shumagin Islands, Bendel Island     <NA>  
## 5    6639 1568_2006-07-31  55.1 -160. Shumagin Islands, Spectacle Island  <NA>  
## 6    6640 1569_2006-07-31  55.1 -160. Shumagin Islands, Nagai Island, Mi… <NA>  
## 7    6641 1569_2006-07-31  55.1 -160. Shumagin Islands, Nagai Island, Mi… <NA>

# Let's address those using Location,
regions = mutate(regions,
       Region = case_when(!is.na(Region) ~ Region,
                          is.na(Region) & str_detect(Location, "Yunaska") ~ "Aleutians Yunaska",
                          is.na(Region) & str_detect(Location, "Aghiyuk") ~ "Aleutians Aghiyuk",
                          is.na(Region) & str_detect(Location, "Shumagin") ~ "Aleutians Shumagin"))

# Check again (should be 0),
filter(regions, is.na(Region)) %>% nrow()

## [1] 0

Now let’s add Region as a factor to our visits df, and show it on a map. Since we’re manipulating our data objects again, let’s also rename visits so that we can keep track of our changes.

# Call our wrangled df version visits.0 and remove the non-version 'visits' if it's there,
if ("visits" %in% ls()) assign("visits.0", visits)
if ("visits" %in% ls()) rm(visits)

# Make Region a factor and join by VisitID,
visits.1 = regions %>%
  select(VisitID, Region) %>%
  distinct() %>%
  mutate(Region = as.factor(Region)) %>%
  left_join(visits.0, ., by = "VisitID")

# Function for regional palette
palette.region = colorFactor(palette = brewer.pal(n = visits.1$Region %>% n_distinct(),
                                                  name = 'Spectral'),
                             domain = factor(visits.1$Region))

# Map
leaflet(visits.1) %>% 
  addTiles() %>%
  addCircleMarkers(lng = ~Lon, lat = ~Lat,
                   color = ~palette.region(Region),
                   fillColor = ~palette.region(Region),
                   popup = ~paste(Region, VisitID),
                   stroke = FALSE) %>%
  addLegend(position = "topright",
            pal = palette.region, values = ~Region,
            title = "Region",
            opacity = 0.75)

Like we saw before, the majority of our samples are located in the Gulf of Alaska (GOA), namely Southeast (SEAk) and Southcentral (SCAK). If we were to look at differences in catach in the GOA, it’d be cool if we could include the Aleutian samples too. Obviously, we’ll want to address those few instances from Aghiyuk, Shumagin, and Yunaska, and ideally we aggregate them into one or two groups. We also see relatively small sample sizes from Bristol Bay and from East Beaufort, which may or may not be prohibitive in conducting meaningful tests. One other consideration is the couple of samples along the Bering Strait coast which are lumped together with a larger concentration of samples from South Chukchi.

Some of our first analyses should see if the more isolated samples are different from the other samples. If not, then we can aggregate. If yes, then we should report any significant tests.

2.4 Temporal distribution

Next we’ll want to take a look at how our samples distribute throughout various time periods. This will give us a good idea of how well our samples may or may not overlap- particularly, we’ll want to see overlaps in years and months for better comparisons of the catch data.

# First create a df specifically containing at time variables
times = visits.1 %>%
  select(VisitID, Date) %>%
  mutate(Year = year(Date),
         Month = month(Date, label = TRUE),
         Week = week(Date),
         Day = yday(Date))

# Function to generate frequency plots for temporal vars
plot_times = function(x) {
  times %>%
    filter(Year > 1980) %>%
    count(!!sym(x)) %>%
    distinct() %>%
    ggplot(., aes(x = !!sym(x), y = n)) +
    geom_col() +
    labs(x = NULL, y = NULL,
         title = paste("No. samples per", x, sep = " "),
         subtitle = "w/o data from 1976") 
}

# Generate plots
map(names(times)[3:6], plot_times)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

By year, we see a large variability in number of samples, as well as a lone 1976 year which I recognize to be Blackburn’s dataset from Lower Cook Inlet. That 1976 is an interesting dataset that could be useful for a look at catches over time in that specific area, but it’ll need to be excluded from any other comparisons. The lowest number of samples occur in 2010, 2020, and 2021. We may want to pay close attention to these years in case they show up as outliers. It seems that 1996-2000 is a period of heavy sampling, followed by the mid 2010s and somewhat active year in 2006.

By month, we see an obvious uptick during the summer which makes sense given winter conditions along most of Alaska’s coast (light, ice, etc.). March and October are interesting months because of the potential for impacts on fishes during years with warmer shoulder seasons. November thru February will likely be removed for our purposes. Our key months look to be June, July and August, with possible inclusion of April, May, and September.

It’s interesting if we look at the by-week and by-day graphs, we see a dip in sampling around mid-June/early-July. Not sure why this is- maybe related to more targeting of juvenile salmonids in early June, or maybe so reduced effort around the 4th of July holiday.

2.5 Year & Month by Region

Let’s now see how our samples overlap by Region and Year, and by Region and Month. For these purposes, let’s combine the Aleutians into one group for now, and remove the 1976 data.

Interesting that each region’s samples center on a different year. At a glance, I see the most overlap in 2006 where the Aleutians and both GOA regions look to have similar sample densities, plus there is a smattering of Chukchi/Beaufort sampling too. In 2009, the Bristol Bay samples may be compared with both GoA regions and the Chukchi/Beaufort, although the samples look more sparse across the regions. And in 2013, the Chukchi/Beaufort has good overlap with both GOA regions again. I see nice periods of comparison between SEAK and SCAK, but there is also a concerning difference in yearly distribution of samples between the two regions. Depending on how important year effects appear- we may need to trim the data for better comparisons.

By month, we see that August contains samples from all regions. June is the center for SEAK and the Aleutians, July for SCAK, and August for Chukchi/Beaufort. Best comparisons from this monthly view may be an August (or lumped July + August) comparison of all Regions. Also, a June (or lumped June thru August) comparison of GOA and the Aleutians would be good. And, we should be able to nicely compare SEAK and SCAK from May to September, potentially including April.

For now, we’ll hold off on excluding any of the data. It’ll be better to revisit these considerations after deciding which research questions to tackle, so subsets of the data can be made specific to each one.

2.6 Mesh Size & Replicates

Let’s make similar visuals as before, but now focus on Mesh Size and Replicates in space and time. First, we can re-map the samples to color by Mesh Size and size of points relative Replicates (large dots = more replicates).

I have an issue discerning among the seven different sizes since those warmer colors do not diverge well. Although there are some patterns to make note of: the largest mesh size of 12.7 only occurs in a small area of SEAK, the 10 mm samples similarly only occur with the 12.7 mm samples (likely the Lundstrom et al. 2022 gear comparison), and the smallest size at 3 mm only occurs in Cook Inlet. I don’t see any obvious spatial patterns in Replicates, although this may not be the best view of it.

Let’s see the frequency of samples for both variables:

# Table view
table(visits.1$MeshSize)

## 
##    3  3.2    6  6.4  9.5   10 12.7 
##  436  821  286   35   82   59  148

table(visits.1$Replicates)

## 
##   1   2   3   4   5   6   7   8   9  10  12 
## 665 533 334 142  60  90  19  21   1   1   1

There are more 12.7 mm samples than I would have guessed which is a shame because if they are completely different then we lose a lot of samples. I am wondering if sizes are accurate- I would think that 3 mm and 3.2 mm are effectively similar in practice. I am also wondering if some projects report, say 3 mm instead of 3.2 mm, because of rounding. I would assume that the nets used were manufactured in the US, which would mean that specifications are given in standard not metric. Let’s do a quick check on how typical mesh sizes convert from inches to millimeters.

# Function based on inch to millimeter conversion
in_to_mm = function (x) {
  return(y = 25.4 * x)
}

# Typical mesh sizes in inches
sizes = c((1/8), (1/4), (3/8), (1/2), (5/8), (3/4))

# Convert to millimeters
in_to_mm(sizes)

## [1]  3.175  6.350  9.525 12.700 15.875 19.050

There is obviously some funny rounding happening in the data. I do not see an issue by lumping mesh sizes that are within one millimeter. Even if the nets were actually different by fractions of a millimeter, I don’t think it would affect the selectivity of fishes all that much- especially considering that most of these researchers were capturing fish to be measured to the nearest millimeter.

Let’s see that map again with this in mind:

Much easier to read! I see some concerning patterns we’ll want to investigate further as we go. The majority of samples appear to be using the 3.2 mm size, including both Arctic groups that seem to use it exclusively. There is a large density of 6.4 mm samples in Northern SEAK and throughout the Aleutians. The same issues are visible concerning the two larger mesh sizes. We should do a comparison of 3.2 mm and 6.4 mm to see if there are any significant differences in catch.

With the simpler color scheme, I am also noticing the difference in Replicates a little easier. Just looking at GOA and Aleutian samples, I see that many of the 6.4 mm samples are also low-Replicate samples. Although, many of our samples only contain one replicate so this may just be something I’m seeing and not actually a concern.

Let’s take a look at how these variables in graphical format. I’ll go ahead and aggregate similar mesh sizes and re-classify the both variables as factors, and add in variables for year and month.

visits.2 = mutate(visits.1,
                  # Replicates as an ordered factor
                  Replicates = as.ordered(Replicates),
                  # Combine similar mesh sizes to the tenth decimal
                  MeshSize = case_when(MeshSize < 4 ~ 3.2,
                                       MeshSize > 4 & MeshSize < 7 ~ 6.4,
                                       MeshSize > 7 & MeshSize < 11 ~ 9.5,
                                       MeshSize == 12.7 ~ 12.7) %>%
                    as.ordered()) %>%
  # Add temporal variables
  left_join(select(times, VisitID, Year, Month, Week, Day), by = "VisitID") %>%
  # Week as an ordered factor (month is already an ordered factor)
  mutate(Week = as.ordered(Week))

Our samples are dominated by 3.2 mm mesh size. It looks to have somewhat even distribution among years, but seems to mimic the density of sampling in SEAK and SCAK in the late 90’s, mid 00’s, and mid 10’s. The 6.4 mm size is much more sporadic with many years missing samples. We also see that larger mesh sizes 9.5 mm and 12.7 mm only occur in a handful of years starting 2013. By month, we see similar distributions among the mesh sizes, with most samples occurring in the summer (Jun - Aug) and shoulder months (Apr, May, Sep).

I do not see the yearly distribution of samples differing all that much among replicates. Except that the 6 - 8 replicates tend to occur later in the dataset. I also do not see a difference in how replicates distirbute among months- similar to what we’ve seen that most samples occur in summer and shoulder months.

Last, let’s see a graph of Mesh Size vs Replicates:

Nothing new here really.

2.7 Preparing data for next steps

2.7.1 Spatial aspects

Spatial effects are probably the most unknown to me and could use deeper investigation. We have obvious regional grouping of samples, with a large majority occurring in SEAK and SCAK. There are interesting pockets of data from the Arctic, Bering Sea, and the Aleutians, but how comparable these “regions” are may be complicated by low sample sizes and confounding factors (Mesh Size or Replicates). I think tackling the GoA is the most useful thing to do at the moment since the other areas are much more isolated and densely sampled.

Questions:

• Basically: how should we treat samples from SEAK, SCAK, and Aleutians?
• Are there regional effects in the GoA?
• Do these effects align with current region labels, or is there a more appropriate break point?
• Can the Aleutian samples be grouped with the GoA, or is it much too different?
• Are the various Aleutian samples appropriate to aggregate?
• Similarly, are there appropriate sub-regional groups that can be defined within SEAK and SCAK?
• What sub-regional groupings make sense? Which, if any, are informative?

After these questions are answered, we can recommend how to treat spatial aspects in the data.

2.7.2 Temporal aspects

Unless we want to do some direct comparisons between 1976 and other years of Kachemak Bay, it just seems too far isolated from the rest of the data to include it in analyses. I will try to address seasonal effects using Day of Year, but move to weeks then months if necessary. We should also keep in mind that April/May and September represent “shoulder” months that are much less sampled compared to June - August period. Some regions may not have samples from shoulder seasons, so reducing the dataset may be required depending on where the comparisons are being made.

Questions:

• Which of Day, Week, Month is most informative as seasonality?

Recommendations:

• Remove 1976 data
• Remove months outside of April - September
• Year as random effect
• Seasonality as fixed effect

2.7.3 Mesh Size

Our 3.2 mm samples are the most represented mesh size occurring in all years of data. The 6.4 mm samples are the next most frequent but occur sporadically from 1999-2003, 2006-2008, 2011-2013, and 2017. The other two sizes are the least frequent and occur in later years of the dataset starting in 2013. By month, we see a similar trend as other variables where most samples occur in Jun - Aug no matter the mesh size.

Questions:

• Is there a significant gear effect among mesh sizes? In particular looking at 3.2 mm vs 6.4 mm

Recommendations

• Drop 9.5 mm and 10.7 mm sizes (likely confounded by isolated regional effects)
• Mesh size as fixed effect

2.7.4 Replicates

Something that’s been obvious is the high density of samples with 1 or 2 replicates. I would expect a larger variability in catch response between one/two replicate samples vs three or more replicate samples. There doesn’t seem to be much of a difference in yearly spread of replicates, except when looking at sample with five or more replicates. Similarly by month, we find most of our samples to occur in Jun - Aug no matter the number of replicates.

Questions

• Can we estimate the number of replicates that results in an informative sample? IOW, is there an ideal number of replicates (e.g., 3) beyond which we do not gain any more meaningful information?

Recommendations

• Number of replicates as fixed effect
• Standardize samples by number of replicates, i.e., average catch per visit
• Standardize samples by another estimate, e.g., sample completeness or species accumulation (if needed)
• Weight samples by how much information on average each level of Replicates provides (if needed)

2.8 Fixing data objects and clean up

Last thing we’ll is fix our visits data based on our recommendations after EDA. We’ll save an original version of the data in case we want it for future analyses.

# Create a new visits object containing a subset of the full data
visits = visits.2 %>%
  # Filter temporal info: remove year 1976 and months outside of April - September
  filter(Year > 1980) %>%
  filter(Month > "Mar" & Month < "Oct") %>%
  # Filter mesh size: remove 9.5 mm and 10.7 mm
  filter(MeshSize < '9.5') %>%
  # Re-order filtered factors to match reduced levels
  mutate(Month = month(Date) %>% as_factor(),
         Week = week(Date) %>% as.ordered(),
         MeshSize = MeshSize %>% as.character() %>% as.ordered())

# Keep a full visits data object as another object
visits.all = visits.2

# Remove unwanted objects in environment
rm(list = ls()[!ls() %in% c('events', 'catch', 'visits', 'visits.all')])

# Reset directory
wd = here()
dirs = wd %>% list.files() %>% str_subset(pattern = "^README|^LICENSE|.md$|.Rproj$", negate = TRUE)
for (i in seq_along(dirs)) {
  name = str_replace_all(dirs[i], "^", "dir.")
  path = str_replace_all(dirs[i], "^", str_c(wd, "/"))
  assign(name, path)
  rm(path, i)
}

3 Exploring catch data

3.1 Summary of catch

Since we’ve reduced our samples based on earlier EDA, we need to update our catch samples to reflect this. This also requires subsetting the events-level data object.

Similar to our visits object, we’ll create a working version of our catch data but keep the version that contains all info.

Now we have a subsetted catch data, let’s take a look at what’s there:

skimr::skim(catch.0)

Data summary

Name	catch.0
Number of rows	116565
Number of columns	10
_______________________
Column type frequency:
character	7
numeric	3
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
Sp_CommonName	0	1.00	7	30	0	185	0
Sp_ScientificName	0	1.00	7	33	0	144	0
Fam_CommonName	0	1.00	4	18	0	33	0
Fam_ScientificName	0	1.00	7	17	0	33	0
LengthType	72266	0.38	2	6	0	5	0
LifeStage	47596	0.59	3	8	0	6	0
Gen_ScientificName	0	1.00	5	15	0	112	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
EventID	0	1.00	4932.27	3013.00	1	3023	4015	7377	10340	▃▇▁▅▃
Length_mm	10322	0.91	87.30	66.48	4	44	70	107	930	▇▁▁▁▁
Count	0	1.00	18.07	602.14	1	1	1	1	99997	▇▁▁▁▁

Looking at our character data, we should go ahead and remove common names simply to reduce clutter. LengthType is not very useful and is missing quite a few observations- we’ll drop this too. We knew LifeStage was an issue when we did our wrangling steps. We can keep it for now, but likely won’t use it til way later. Hopefully, we can infer LifeStage down the road by using the length information.

Let’s add our VisitID variable so we analyse our catch at the sample level (instead of replicate).

catch.1 = catch.0 %>%
  # Add VisitIDs
  left_join(select(events, VisitID, EventID), by = "EventID") %>%
  # Remove unwanted parameters
  select(-c(Sp_CommonName, Fam_CommonName, LengthType)) %>%
  # Tidy parameters
  select(VisitID, EventID,
         Fam_ScientificName, Gen_ScientificName, Sp_ScientificName,
         Count, Length = Length_mm, LifeStage)

3.2 Calculate metrics of catch

3.2.1 Occurrence and rare taxa

Previously in the wrangle step, we defined rare as occurring in less than 0.25% of samples. However, we may start with a more conservative approach to see if removing more rare species is necessary. It makes sense to me to remove species with only 1 or 2 occurrences for now. This is equivalent to removing species occurring in less than 0.15% of samples.

# Create a tibble for species occurrences
sp.occur = catch.1 %>% 
  # For each species within each visit, Presence = 1
  summarise(Presence = n_distinct(Sp_ScientificName), .by = c(VisitID, Sp_ScientificName)) %>%
  # For each species, total the number of presences
  summarise(Occurrence = sum(Presence), .by = Sp_ScientificName) %>%
  # Calculate occurrence as a percentage
  mutate(Perc_Occurrence = round(Occurrence / n_distinct(events$VisitID) * 100, 2)) %>%
  # Add rare category (subject to change)
  mutate(Rare = ifelse(Occurrence < 3, "Yes", "No"))

# Create a tibble for total presences per sample
n.occur = catch.1 %>% 
  # For each species within each visit, Presence = 1
  summarise(Presence = n_distinct(Sp_ScientificName), .by = c(VisitID, Sp_ScientificName)) %>%
  # For each visit, total the number of presences
  summarise(n_sp = sum(Presence), .by = VisitID)

# Create a vector of the rare species
rare.species = filter(sp.occur, Rare == "Yes") %>% pull(Sp_ScientificName)

# Remove rare taxa and update the catch tibble
catch.2 = filter(catch.1, !Sp_ScientificName %in% rare.species)

3.2.2 Average abundance per visit

Because our samples (visits) contain a variable number of replicates (events), I decide to use an averaging method to calculate abundance. Here, we use a simple mean (i.e., total count / number of replicates) then I round to the nearest upper integer. Note that this rounding method gives more representation of rarer taxa. For example, say we have a sample with three replicates and caught species X at 0, 0, and 1 frequency. We get an average of 0.33 which results in an abundance of 1 after rounding up. In another scenario let’s say we catch species Y in the same sample at 20, 22, and 400 frequency. We have an average of 147.3 resulting in an abundance of 148 after rounding up. This is our first instance of data standardization, but we may want to further standardize and/or transform abundance data depending on EDA.

# Create a tibble for avg sp abundance per visit
abun = catch.2 %>%
  select(VisitID, EventID, Sp_ScientificName, Count) %>%
  summarise(Abundance = (sum(Count) / n_distinct(EventID)) %>% ceiling(), .by = c(VisitID, Sp_ScientificName))

skimr::skim(select(abun, Abundance))

Data summary

Name	select(abun, Abundance)
Number of rows	13031
Number of columns	1
_______________________
Column type frequency:
numeric	1
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Abundance	0	1	76.21	833.45	1	1	3	10	45557	▇▁▁▁▁

3.2.3 Size information

# Create a tibble for containing length data per species and visit
size = select(catch.2, VisitID, EventID, Sp_ScientificName, Length)

3.3 Diversity

3.3.1 Choosing indices

# Create a tibble of diversity indices per sample
diversity = abun %>%
  group_by(VisitID) %>%
  summarise(S = specnumber(Sp_ScientificName), # Richness
            H = diversity(Abundance, index = "shannon"), # Shannon
            # q1 = exp(H), # exp(Shannon)
            D = diversity(Abundance, index = "simpson"), # Simpson
            # q2 = diversity(Abundance, index = "invsimpson"), # Inverse-Simpson
            J = diversity(Abundance, index = "shannon") / log(S)) # Pielou's evenness

3.3.2 Plot diversity across variables

# Function to plot diversity indices
plot_diversity = function(x) {
  dat = visits %>%
    select(VisitID, !!sym(x)) %>%
    left_join(diversity, by = "VisitID")
  plot_indices = function(dat) {
    plot_list = list()
    D.indices = names(dat)[3:6]
    for (i in 1:length(D.indices)) {
      pluck(plot_list, D.indices[i]) = ggplot(dat) +
                                        aes(x = !!sym(x), y = !!sym(D.indices[i])) +
                                        geom_jitter(alpha = 0.5)
    }
    out = ggarrange(plotlist = plot_list,
                    ncol = 2, nrow = 2,
                    common.legend = TRUE)
    return(out)
  }
  plot_indices(dat) %>% return()
}

# Define which variables to plot diversity against and save as a vector 
vars = c("Replicates", "MeshSize", "Lon", "Year", "Month")

# For each defined variable, arrange plots of diversity
map(vars, suppressMessages(plot_diversity))

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

3.3.3 Differences in richness

# Visualize richness by replicates
left_join(visits, diversity, by = "VisitID") %>%
  boxplot(S ~ Replicates, data = .)

# remove replicates = 9, 10, 12 bc of single obs
tmp = left_join(visits, diversity, by = "VisitID") %>%
  mutate(Replicates = as.numeric(Replicates)) %>%
  filter(Replicates < 9)

# analysis of variance
aov.S.reps = aov(S ~ as.factor(Replicates), data = tmp)
summary(aov.S.reps)

##                         Df Sum Sq Mean Sq F value Pr(>F)    
## as.factor(Replicates)    7  14067  2009.5   80.91 <2e-16 ***
## Residuals             1469  36485    24.8                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# pairwise test
TukeyHSD(aov.S.reps) %>% broom::tidy() %>% select(-term) %>% knitr::kable(digits = 3)

contrast	null.value	estimate	conf.low	conf.high	adj.p.value
2-1	0	3.133	2.184	4.081	0.000
3-1	0	7.451	6.241	8.661	0.000
4-1	0	6.231	4.732	7.729	0.000
5-1	0	8.400	5.474	11.325	0.000
6-1	0	9.442	7.079	11.805	0.000
7-1	0	4.995	0.584	9.406	0.014
8-1	0	8.434	4.908	11.959	0.000
3-2	0	4.318	3.056	5.581	0.000
4-2	0	3.098	1.557	4.639	0.000
5-2	0	5.267	2.320	8.214	0.000
6-2	0	6.309	3.919	8.700	0.000
7-2	0	1.862	-2.563	6.287	0.907
8-2	0	5.301	1.757	8.844	0.000
4-3	0	-1.220	-2.935	0.494	0.377
5-3	0	0.949	-2.093	3.990	0.981
6-3	0	1.991	-0.515	4.497	0.236
7-3	0	-2.456	-6.945	2.033	0.713
8-3	0	0.983	-2.640	4.605	0.992
5-4	0	2.169	-0.998	5.336	0.429
6-4	0	3.211	0.554	5.868	0.006
7-4	0	-1.236	-5.811	3.339	0.992
8-4	0	2.203	-1.526	5.932	0.625
6-5	0	1.042	-2.615	4.699	0.989
7-5	0	-3.405	-8.624	1.814	0.496
8-5	0	0.034	-4.462	4.530	1.000
7-6	0	-4.447	-9.373	0.479	0.112
8-6	0	-1.008	-5.161	3.144	0.996
8-7	0	3.439	-2.139	9.016	0.571

# residuals
residuals = residuals(aov.S.reps)
plot(residuals)

par(mfrow = c(2, 2))
plot(aov.S.reps)

par(mfrow = c(1, 1))
# check for homogeneity
shapiro.test(residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals
## W = 0.95897, p-value < 2.2e-16

# result: violation of homogeneity

# kruskal-wallis test
kruskal.test(S ~ Replicates, data = tmp)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  S by Replicates
## Kruskal-Wallis chi-squared = 422.05, df = 7, p-value < 2.2e-16

# one way t-test not assuming equal variances
oneway.test(S ~ Replicates, data = tmp)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  S and Replicates
## F = 82.534, num df = 7.000, denom df = 97.828, p-value < 2.2e-16

# pairwise test
results = pairwise.t.test(tmp$S, tmp$Replicates,
                          p.adjust.method = "BH",
                          pool.sd = FALSE)
results

## 
##  Pairwise comparisons using t tests with non-pooled SD 
## 
## data:  tmp$S and tmp$Replicates 
## 
##   1       2       3       4       5       6       7      
## 2 < 2e-16 -       -       -       -       -       -      
## 3 < 2e-16 < 2e-16 -       -       -       -       -      
## 4 < 2e-16 2.9e-06 0.11531 -       -       -       -      
## 5 3.6e-08 7.5e-05 0.45956 0.11531 -       -       -      
## 6 < 2e-16 1.7e-11 0.02676 0.00083 0.45956 -       -      
## 7 0.01418 0.30978 0.18341 0.46741 0.11531 0.02886 -      
## 8 1.3e-06 0.00035 0.45956 0.11742 0.98221 0.46741 0.11531
## 
## P value adjustment method: BH

# plot t.test results
ggplot(tmp) +
  aes(x = Replicates, y = S, group = Replicates) +
  geom_boxplot() +
  geom_signif(test = "t.test",
              comparisons = list(c(1, 2),
                                 c(2, 3),
                                 c(2, 4),
                                 c(3, 4),
                                 c(3, 5)),
              y_position = c(33, 36, 39, 43, 46),
              map_signif_level = TRUE)

3.4 Species EDA plots

Walleye Pollock

Starry Flounder

Great Sculpin

Pacific Herring

Coho Salmon

Dolly Varden

3.5 MV Ordinations

3.5.1 Create catch matrix based on average abundances

# Create a data.frame in sample x species format (e.g., vegan)
abun.mat = abun %>%
  pivot_wider(names_from = Sp_ScientificName,
              values_from = Abundance,
              values_fill = 0) %>%
  column_to_rownames(var = "VisitID")

3.5.2 Standardize catch

# presence/absence
pa = decostand(abun.mat, method = "pa")

# 4th root
rt4 = mutate(abun.mat, across(everything(), ~ .x^(1/4)))

# robust CLR
rclr = decostand(abun.mat, method = "rclr")

3.5.3 Calculate distances

# Jaccard distances on P/A data
jac = vegdist(as.matrix(pa), method = "jaccard")

# Bray-Curtis distance on 4th rt transformed data
bc = vegdist(as.matrix(rt4), method = "bray")

# Euclidean distance on RCLR transformed data, i.e., robust aitchison distance
rait = vegdist(as.matrix(rclr), method = "euclidean")

3.5.4 Ordinations

nmds.jac = metaMDS(jac)
plot(nmds.jac, type = "t")

nmds.bc = metaMDS(bc)
plot(nmds.bc, type = "t")

## plot using rda() abd biplot()
# pca.rclr = rda(as.matrix(rclr))
# biplot(pca.rclr, display = c("sites", "species"), type = c("text", "points"))

mds.rclr = metaMDS(rait)
plot(mds.rclr, type = "t")

3.5.5 Look at the outlier sample

filter(visits, VisitID == "846_2015-08-15") %>%
  left_join(catch.2, by = "VisitID") %>%
  glimpse()

## Rows: 1
## Columns: 19
## $ VisitID            <chr> "846_2015-08-15"
## $ Replicates         <ord> 1
## $ Date               <date> 2015-08-15
## $ Lat                <dbl> 71.23332
## $ Lon                <dbl> -155.7281
## $ MeshSize           <ord> 3.2
## $ Cluster            <fct> 101
## $ Region             <fct> Chukchi/Beaufort
## $ Year               <dbl> 2015
## $ Month              <fct> 8
## $ Week               <ord> 33
## $ Day                <dbl> 227
## $ EventID            <dbl> 5396
## $ Fam_ScientificName <chr> "Coregoninae"
## $ Gen_ScientificName <chr> "Coregonus"
## $ Sp_ScientificName  <chr> "Coregonus laurettae"
## $ Count              <dbl> 2
## $ Length             <dbl> NA
## $ LifeStage          <chr> NA